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Abstract In the Heston stochastic volatility model, the transition probability of the 
variance process can be represented by a non-central chi-square density. We focus on the 
case when the number of degrees of freedom is small and the zero boundary is attracting 
and attainable, typical in foreign exchange markets. We prove a new representation for 
this density based on sums of powers of generalized Gaussian random variables. Further 
we prove Marsagiia's polar method extends to this distribution, providing an exact 
method for generalized Gaussian sampling. The advantages are that for the mean- 
reverting square-root process in the Heston model and Cox-Ingersoll-Ross model, we 
can generate samples from the true transition density simply, efficiently and robustly. 
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1 Introduction 

A popular choice to model the evolution of an asset price and its stochastic volatility 
is the correlated Heston model [23] : 

dS t = pSt At+sJVtSt {pdWt 1 + a/1-P 2 dW?), 
dV t = k{9 ~ V t ) dt + e^VtdWl. 

Here (W 1 ,!^ 2 ) is a standard two-dimensional Wiener process, the parameters p, n, 
9 and e are positive constants, and p £ (— 1, 1). The process S represents the price 
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process of an underlying asset, for example a stock index or an exchange rate, and V is 
the variance process of the log returns of S, modelled as a mean-reverting square-root 
process, the well-known Cox-Ingersoll-Ross process. 

The transition probability density of the variance process V is known explicitly, it 
can be represented by a non-central chi-square density. Depending on the number of 
degrees of freedom v := AkO/e 2 , there are fundamental differences in the behaviour of 
the variance process. If v is larger or equal to 2, the zero boundary is unattainable; if it 
is smaller than 2, the zero boundary is attracting and attainable. At the zero boundary 
though, the solution is immediately reflected into the positive domain. 

A number of successful simulation schemes have been developed for the non- 
attainable zero boundary case. There are schemes based on implicit time-stepping 
integrators, see for example Alfonsi [3][4] and Kahl and Schurz 33 . Other time dis- 
cretization approaches involve splitting the drift and diffusion vector fields and eval- 
uating their separate flows (sometimes exactly) before they are recomposed together, 
typically using the Strang ansatz. See for example Higham, Mao and Stuart 26 , Ni- 
nomiya and Victoir [45] and Alfonsi [I]. However, the splitting methods and the implicit 
methods so far only apply in the non-attracting zero boundary case. 

Other direct discretization approaches, that can be applied to the attainable and 
unattainable zero boundary case are based on forced Euler-Maruyama approxima- 
tions and typically involve negativity truncations; some of these methods are positivity 
preserving. See for example Deelstra and Delbaen |14| . Bossy and Diop [5] and also 
Berkaoui, Bossy and Diop [7], Lord, Koekkoek and Van Dijk [39], as well as Higham and 
Mao 25 , among others. These methods all converge to the exact solution, but their rate 
of strong convergence is difficult to establish. However, this said, the leading method 
in this class — the full truncation method of Lord, Koekkoek and Van Dijk [3U] — has in 
practice proved highly effective. 

Exact simulation methods typically sample from the known non-central chi-square 
distribution xiW f° r the transition probability of the variance process V (see Cox, 
Ingersoll and Ross [I3j and Glassermann [191 Section 3.4]). Broadie and Kaya [10] pro- 
posed sampling from xiW as follows. When v > 1, xiW = (N(0, ^A)) + Xv-li so 
such a sample can be generated by a standard Normal sample and a central chi-square 
sample. When < v < 1, such a sample can be generated by sampling from a Poisson 
distribution with mean A/2, and then sampling from a central X2N+v distribution. To 
simulate the asset price Broadie and Kaya integrated the volatility process to obtain an 
expression for J y/VrdWr- They substituted that expression into the stochastic differ- 
ential equation for In St- The most difficult task left is then to simulate J V T dr on the 
global interval of integration conditioned on the endpoint values of V; see Smith [53j . 
The Laplace transform of the transition density for this integral is known from results 
in Pitman and Yor 47 . Broadie and Kaya used Fourier inversion techniques to sample 
from this transition density. Glasserman and Kim [20] on the other hand, showed that 
linear combinations of series of particular gamma random variables exactly sample this 
density. They used truncations of those series to generate suitably accurate sample ap- 
proximations. Anderson 5 suggested two approximations that make simulation of the 
Heston model very efficient. The first was, after discretizing the time interval of integra- 
tion for the price process, to approximate / Vr dr on the integration subinterval by a 
simple quadrature. This would thus require non-central xt(X) samples for the volatility 
at each timestep. Hence the second, was to approximate and thus efficiently sample, 
the X?(A) distribution — in two different ways depending on the size of A. Haastrecht 
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and Pelsser [22] have recently introduced a rival xiW sampling method to Andersen's. 
Moro and Schurz [44] have also successfully combined exponential splitting with exact 
simulation. Dyrting [K] outlines and compares several different series and asymptotic 
approximations for non-central chi-square distribution. 

There are also numerous approximation methods based on the corresponding Fokker- 
Planck partial differential equation. These can take the form of Fourier transform 
methods — see Carr and Madan [1 1] . Kahl and Jackel [32] or Fang and Oosterlee |16l 
I17| for example — or some involve direct discretization of the Fokker-Planck equation. 

We focus on the challenge of the attainable zero boundary case and in particular on 
the case when v <C 1, typical of FX markets and long-dated interest rate markets as 
remarked in Andersen [5] . The method we propose follows the lead of Andersen [5] . It is 
based on efficiently simulating the known non-central xiW transition density for the 
volatility required for each timestep of Andersen's integration method for the Heston 
model. More precisely, we suggest an exact simulation method for the non-central XuW 
density as follows. A non-central xt{X) random variable can be generated from a central 
X2N+1/ random variable with N chosen from a Poisson distribution with mean A/2. 
Further, a X2N+v random variable can be generated from the sum of squares of 2N 
independent standard Normal random variables and an independent central xi random 
variable. So the question we now face is how can we efficiently simulate a central xi 
random variable, especially for v < 1? Suppose that v is rational and expressed in 
the form v = p/q with p and q natural numbers. We show that a central Xv random 
variable can be generated from the sum of the 2gth power of p independent random 
variables chosen from a generalized Gaussian distribution N(0, 1, 2q), where a N(0, 1, q) 
distribution has density 

/N(0,i,9) (*) : = 2 l+l/ g 9 r(1/g) ■ ^pHW)' 

where x € M and F(-) is the standard gamma function. How can we sample from 
a N(0, l,2g) distribution? Our answer lies in generalizing Marsaglia's polar method 
for pairs of independent standard Normal random variables. Indeed we generate 2q 
uniform random variables U = (U\, . . . , U2 q ) over [—1, 1], and condition on their 2qih 
norm ||[/||2g, being less than unity. Then we prove that the 2q random variables U ■ 
(— 21og H^l^g) 2 VII^1l2g are independent N(0, 1, 2q) random variables. We provide a 
thorough comparison, of our generalized Marsaglia polar method for sampling from 
the central xi distribution, to the acceptance-rejection method of Ahrens and Dieter 
(see Ahrens and Dieter [2] and Glasserman [19] 1. 

The Cox-Ingersoll-Ross process, which has a non-central chi-squared transition 
probability, can thus be exactly simulated by the approach just described, which we 
will call the Marsaglia generalized Gaussian method (MAGG). The advantages of this 
approach are that for the mean-reverting variance process in the Heston model, we can 
generate high quality samples simply and robustly. The disadvantage of MAGG when 
simulating the Cox-Ingersoll-Ross process is that the degrees of freedom v must be 
rational (however, in practice, this is typically fulfilled). We demonstrate our method 
in the computation of option prices for parameter cases that are considered in Ander- 
sen [5] and Glasserman and Kim [2(J| and described there as challenging and practically 
relevant. To summarize, we: 

— Prove that a central chi-squared random variable with less than one degree of free- 
dom, can be written as a sum of powers of generalized Gaussian random variables; 
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— Prove a new method — the generalized Marsaglia polar method — for generating gen- 
eralized Gaussian samples; 

— Establish a new simple, exact, unbiased and efficient method for simulating the 
Cox-Ingersoll-Ross process, for an attracting and attainable zero boundary, and 
thus establish a new simple method for simulating the Heston model. 

Our paper is organised as follows. In Section [2] we present our new method for 
sampling from the non-central chi-squared distribution based on sampling from the 
generalized Gaussian distribution. We include a thorough numerical comparison with 
the acceptance-rejection method of Ahrens and Dieter [2j. We apply our MAGG method 
to the Heston model in Section[3] We compare its accuracy and efficiency to the method 
of Andersen [5]. Finally in Section [4] we present some concluding remarks. 

2 Non-central chi-square sampling 

We present our new theoretical probabilistic results. We begin by introducing the 
generalized Gaussian distribution and show that central xi random variables can be 
represented by sums of powers of generalized Gaussian random variables. To generate 
central xi samples in this way, we require an efficient method for generating gener- 
alized Gaussian samples. We provide such a method in the form of a generalization 
of Marsaglia's polar method for standard Normal random variables. Then, to put our 
approach on a firm practical footing, we compare its efficiency to the most well-known 
exact xi sampling method, the acceptance-rejection method of Ahrens and Dieter [2J. 
Rounding off this section, we return to our original goal, and explicitly state our general 
algorithm for non-central xiW- This will be the algorithm we use in our application 
to the Heston model in Section [3] 

2.1 Generalized Gaussian distribution 

We prove that random variables with a central Xv distribution, especially for v < 1, 
can be represented by random variables with a generalized Gaussian distribution. 

Definition 1 (Generalized Gaussian distribution) A generalized N(0, l,q) ran- 
dom variable, for q ^ 1, has density 



where igl and r(-) is the standard gamma function. 

See Gupta and Song |21j . Song and Gupta |54] , Sinz, Gerwinn and Bethge |52] and 
Sinz and Bethge |51] for more details on this distribution and its properties. 

Theorem 1 (Central chi-square from generalized Gaussians) Suppose Xi ~ 
N(0, 1,2^) are independent identically distributed random variables for i — l,...,p, 
where q ^ 1 and p £ N. Then we have 




p 
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Proof If X ~ N(0, 1, 2q), then we see that 

r\x\ 1/2q 

where £ = M ■ Hence we deduce that \X\ 2q ~ Xi/ g - Now using that the sum of 

p independent identically distributed Xi/q random variables have a Xp/ q distribution 
establishes the result. □ 



2.2 Generalized Marsagiia polar method 

If we intend to use N(0, 1, 2q) samples to generate Xp/ q samples, we need an accurate 
and efficient method for sampling from a generalized Gaussian distribution. To this 
end we generalize Marsagiia's polar method for pairs of independent standard Normal 
random variables (see Marsagiia [421). 

Theorem 2 (Generalized Marsagiia polar method) Suppose for some q £ N 

that U\, , Uq are independent identically distributed uniform random variables over 

[—1,1]. Condition this sample set to satisfy the requirement \\U\\ q < 1, where \\U\\ q 
is the q-norm of U = (U\, ■ ■ ■ , U q ). Then the q random variables generated by U ■ 
(— 2 log ||C^||q) 1//9 /||!7|| g are independent N(0,l,q) distributed random variables. 

Proof Suppose for some q £ N that U = (U\,...,U q ) are independent identically 
distributed uniform random variables over [—1,1], conditioned on the requirement that 
|| U || q < 1. Then the scalar variable 

Z := (-21og ire 1/9 >0 

is well defined. Let / denote the probability density function of U given \\U\\ q < 1; it is 
defined on the interior of the g-sphere, §<j(l), whose bounding surface is \\U\\ q = 1. We 
define a new set of q random variables W = (Wj., ■ • ■ , W q ) by the map G : §q(l) — > R p 
where G: U h-> W is given by 



Go U = 




Note that the inverse map G 1 : W — > S q (l) is well defined and given by 

where we note that in fact Z = ||W||g which comes from taking the g-norm on each 
side of the relation W = G(U). 

We wish to determine the probability density function of W . Note that if Q C R* 3 , 

P(W e ri) = P(u e G~\n)) 

— f o udu 

(/ o G _1 o w) ■ |det(DG~ 1 o w) I dw, 
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where for w = . . . ,w p ) £ O, the quantity DG 1 o w denotes the Jacobian trans- 
formation matrix of G _1 . Hence the probability density function of W is given by 

(foG^ 1 ow) ■ dct(DG _1 ow)\. 

The Jacobian matrix and its determinant are established by direct computation. For 
each i, k = 1, . . . , q we see that if we define g(z) := —(1/2 + 1/ z q ), then 



dGj/ = exp(-zV2g) {/ . x ^ nA Ufi<|9 _^ 



(<*ifc + 3(2) ■ (sgn(w l ) • |^| 9 1 ) ■ w fe J , 
where 5n- is the Kronecker delta function. If we set 

v = (sgn(uii) ■ ItoiI 9-1 , . . . ,sgn(w g ) ■ |ui g | 9_1 ) T 

then we see that our last expression generates the following relation for the Jacobian 
matrix (here I q denotes the 5x9 identity matrix): 

(DG -1 on) — Iq + g(z) ■ v ui T . 



exp(~zi/2q) 

From the determinant rule for rank-one updates — see Meyer [431 p. 475] — we see that 
the determinant of the Jacobian matrix is given by 

det(DG~ 1 o w ) = eXp( ~ Z?/2) ■ (1 + g( z ) W T v) 
= exp^V2)_ (l + ff(z) ^ 

= -£exp(-*»/2), 

where we used the definition for g(z) in the last step. Noting that vol(8g(l)J = 2 q ■ 
(r(l/q)) q /q q we have 

(/ oG^ow). |det(DG~ 1 ow)| = • exp(-2 9 /2). 

V ' 1 V ' l 21+ 1 (r(l/q)) q V ' ' 

This is the joint probability density function for q independent identically distributed 
g-generalized Gaussian random variables, establishing the required result. □ 



2.3 Comparison with acceptance-rejection method 

We compare our approach for generating central Xi> samples to the acceptance-rejection 
algorithm of Ahrens and Dieter [2] for the case v < 2. In particular we use the form 
of this acceptance-rejection algorithm outlined in Glasserman [191 pp. 126-7]. The 
acceptance-rejection algorithm is based on a mixture of the prior densities (v/2) a;"/ 2-1 
on [0,1] and exp(l — x) on (l,oo), with weights e/(e + v/2) and (f/2)/(e + v/2), 
respectively; here e = exp(l). This method generates one Xv random variable with 
probability of acceptance 

. {v/2)r(v/2)e 
PAD "= ,,/2 + c ' 
In this method, the number of degrees of freedom v can be any real number. 

For our generalized Marsaglia polar method, we restrict ourselves to case when the 
number of degrees of freedom is rational, i.e. v = p/q with p, q £ N. The algorithm for 
generating central xi samples is as follows. 
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Algorithm 1 (Central chi-square samples) To produce an exact Xp/ q sample: 

1. Generate 2q independent uniform random variables over [—1,1]: U — (U\, . . . , U% q ). 

2. If ||!7||2g < 1 continue, otherwise repeat Step 1. 

3. Compute Z = U ■ (-2 log \\U\\l q ) 1/2q /\\U\\ 2q - This gives 2q independent N(0, 1, 2q) 
distributed random variables Z — (Z\, . . . , Zi q \ 

4. Compute Z 2q -\ h Zp 9 . 

In the second step, the probability of accepting U\, . . . , Ui q is given by the ratio of 
the volumes of S 2 ,(l) and [-1, l] 2q : 



2 q 

Ma 



r(i/2 g ) V 

2q I 



Note for q = 1, the probability of acceptance is 0.7854. Further as q — > oo we have 
-PMar — * ex P( — 7) ~ 0.5615. Here 7 is the Euler-Mascheroni constant and we have used 
that P(z) ~l/z-7asz^ + . 

In practice we will need to generate a large number of samples. For the generalized 
Marsaglia polar method, in each accepted attempt, we generate 2q generalized Gaus- 
sian random variables. Of these, p random variables are used to generate a X^,/ q random 
variable. The number of attempts until the first success has a geometric distribution 
with mean 1/PMar- Hence the expected number of steps to generate 2q/p indepen- 
dent Xp/ q random variables is thus 1/PMar- F° r the acceptance-rejection method, the 
expected number of attempts to generate 2q/p independent Xp/ q distributed random 
variables is {2q/p) ■ (1/Pad)- 

The first natural question is whether the expected number of attempts for the gen- 
eralized Marsaglia polar method is less than that for the acceptance-rejection method. 
In other words, to generate 2q/p random variables, is l/^Mar ^ (2?/p) ■ (1/Pad)? Or 
equivalently, when does p/q ^ 2PMar/^ D AD hold? We examine the right-hand side more 
carefully; set z := l/2q, so < z < 1/2. Then we have 

D Mar _ (_ r.^.A^ v l 2 + c 



(zr(z)) 



Pad V V /; 0V2)i>/2)e' 

Note z and v/2 are independent. A lower bound for ^f(z)) is exp(— 7) ~ 0.5615 for 
< z < 1/2, whilst a lower bound for (v/2 + e)/((f/2) r{v/2) e) is 1 for < v < 2. 
Hence 2PMar/-PAD > 1 an d so for p/q < 1, the expected number of attempts for the 
generalized Marsaglia method is less than that for the acceptance-rejection method. 
We further note that the expected number of attempts for the generalized Marsaglia 
method to generate 2q/p chi-square samples is bounded by its value for q = 1 and 
the limit as q — > 00, more precisely the expected number of attempts is 1/PMar 6 
(1.2732, 1.7811). In contrast in the acceptance-rejection method, the expected number 
of attempts to generate 2q/p samples is (2q/p)(l/ P^), which is unbounded. 

The second natural question is how do the two algorithms perform in practice? The 
issue that immediately surfaces is that the generalized Marsaglia method is restricted to 
rational numbers. However, in practical applications this is not restrictive. In Figure [T] 
the upper panel shows the the CPU time needed by both methods to generate 10 6 
central xi samples for the values u = n ■ 10 _m where n — 1, . . . , 9 and m = 2,3, 4. 
We observe that the Ahrens and Dieter acceptance-rejection method roughly requires 
the same CPU time to generate central xi samples for these values of v. It is slower 
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Fig. 1 CPU time versus the number of degrees of freedom v for the generalized Marsaglia 
method and the Ahrens and Dieter acceptance-rejection method. The upper panel shows the 
CPU time needed to generate 10 samples for the values v = n ■ 10 — m where n = 1, ... ,9 
and m = 2,3, 4. The lower panel shows the CPU time to need to generate 10 4 samples, 
simultaneously for two sets of v values given to three significant figures, namely u = (1 -f Wi) ■ 
10" 4 for m = 0, 1, . . . , 1000 and v = 0.101 + m ■ 10" 3 for m = 0, 1, . . . , 299. 
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than for the generalized Marsaglia method. However the generalized Marsaglia method 
shows more variation in the CPU time required. In particular for example, for values 
of v equal to 3, 6, 7 and 9 times 10 _m for all m = 2, 3, 4, it takes longer to generate 
central Xv samples than for the other v values. This is due to the fact that as rational 
numbers, with denominators as powers of 10, they do not simplify nicely to what 
might be considered the optimal format for sampling with this method, namely 1/q. 
For values of v which cannot be reduced to this optimal format, we need to sum over 
a number of generalized Gaussian samples to produce a central xi sample. However 
any decimal with a finite number of significant figures can be written as the sum of 
fractions of powers of 10. Further a central Xu random variable can be constructed 
by adding central Xu t random variables for which v\ + ■ ■ ■ + — v. Suppose indeed, 
we generate a central Xv sample by adding xti samples where the z/j are fractions of 
powers of 10 that generate each significant figure. From the upper panel in Figure[T] we 
observe that provided v does not have too many significant figures, then on average, 
the generalized Marsaglia method will be more efficient than the acceptance-rejection 
method. This is in fact confirmed in the lower panel in Figure[T] There we show how the 
CPU time varies with the number of degrees of freedom, when v is given to 3 significant 
figures. Even with the number of degrees of freedom given to 4 significant figures, we 
can see from the lower panel in Figure [T] that the generalized Marsaglia method will 
still be more efficient on average. Parameter values such as the number of degrees of 
freedom v are often determined by calibration. Since these are often quoted to only 
2 or 3 significant figures, the generalized Marsaglia method would be the method of 
choice. We return to this issue in our conclusion section. 

All simulations were run in Matlab, whose Profiler feature reveals that for the 
generalized Marsaglia method most CPU time is spent on computing ||?7||2g in Step 2 
and the sum in Step 4. In the acceptance-rejection method most time is spent on the 
decision processes required for choosing which of the mixture of prior densities to use. 



2.4 Non-central chi-squared sampling 

We return to exact simulation of the non-central xS(A) distribution. A xiW random 
variable can be generated as follows. Choose a random variable N from a Poisson dis- 
tribution with mean A/2. Then a sample generated from a central X2N+v distribution 
is in fact a xtW sample. In other words we have 

oo 

See for example Johnson [29] or Broadie and Kaya [10] • Hence we are left with the 
problem of how to sample from a x\n+v distribution. If we generate 2N independent 
standard Normal random variables, say Y\,..., V2JV1 an d an independent Xu random 
variable, say Z, then Yj +■ ■ --\-Y^-\-Z ~ X2N+v Putting all the components together 
we arrive at the following simple algorithm. We assume v = p/q with p and q natural 
numbers. 

Algorithm 2 (Non-central chi-square samples) To produce an exact Xp/ q W 
sample: 
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1. Use Algorithm \T\ to generate 2q independent N(0, l,2g) distributed random vari- 
ables Z=(Z 1 ,...,Z 2q ). 

2. Generate Poisson distributed random variable N with mean A/2. 

3. Generate 2N standard Normal random variables, call them Yi, . . . , l^jy- 

4. Compute Y 2 + ■ ■ ■ + Y$ N + Z 2q + ■ ■ ■ + Z 2 p q ~ X 2 p/q W- 

Note that if p < 2q then we can use the remaining N(0, 1, 2q) random variables we 
generate in Step 3, the next time we need to generate a Xp/ q W sample. In practice 
we don't really need to consider the case p ^ 2q, but for the sake of completeness, we 
would simply generate p — 2q more N(0, 1, 2q) samples by repeating Steps 1-3. 

Remark 1 (Poisson sampling) As Haastrecht and Pelsser ;22] note, the mean /i = A/2 
of the Poisson distribution we wish to sample from is small. We can efficiently draw 
a sample from such a Poisson distribution by inverting its distribution function over 
a uniform random variable — see Knuth [361 p. 137] or Ahrens and Dieter pQ. The 
algorithm is as follows. Calculate exp(— fi). Generate uniform random variables (over 
[0, 1]) say U\, U2, ■ ■ ■ until U\ ■ U2 ■ •■ U m ^ exp(— p). Set N •<— m — 1. On average this 
algorithm requires the generation of p + 1 uniform variates. 



3 Application: the Heston model 

The Heston model (Heston [23]) is a two- factor model, in which the first component 
S describes the evolution of a financial variable such as a stock index or exchange 
rate, and the second component V describes the stochastic variance of its returns. The 
Heston model is given by 

dSt = pSt dt + p^VtSt dWt 1 + \/l-P 2 VVtSt dW t 2 , 
dV t = k(6 ~ V t ) dt + esJVtdWh 

where and W 2 are independent scalar Wiener processes. The parameters /i, k, 6 
and e axe all positive and p G (— 1, 1). In the context of option pricing, a pricing measure 
must be specified. We assume here that the dynamics of S and V as specified above are 
given under the pricing measure. For a discussion and derivation of various equivalent 
martingale measures in the Heston model see for example Hobson [24]. By the Yamada 
condition this model has a unique strong solution. In particular, the volatility V is non- 
negative, and the stock price S, as a pure exponential process, is positive. Without loss 
of generality we suppose /j, — 0. We explain how we simulate the variance process V 
for the case of an attracting and attainable zero boundary first. Then we discuss how 
we approximately simulate the asset price process S. Finally we present some explicit 
simulation results. 



3.1 Variance process simulation 

The variance process Vt, generated by the scalar stochastic differential equation above, 
is modelled as a mean-reverting process with mean 8, rate of convergence k and square 
root diffusion scaled by e. It is known as the Cox-Ingersoll-Ross process (see Cox, 
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Ingersoll and Ross [13] who modelled the short rate of interest using this process). We 
define the degrees of freedom for this process to be 

v := 4k6»/e 2 . 

When the process Vt can be reconstructed from the sum of squares of v Ornstein- 

Uhlenbeck processes; hence the label of degrees of freedom. When v < 2 the zero bound- 
ary is attracting and attainable, while when v ^ 2, the zero boundary is non-attracting. 
These properties are immediate from the Feller boundary criteria, see Feller [18] . These 
are based on inverting the associated stationary elliptic Fokker-Planck operator, with 
boundary conditions, and can be found for example in Karlin and Taylor 35 . 

Here we focus on the challenge of v < 2 and in particular cases when v <C 1 typi- 
cal of FX markets (Andersen 5!). Importantly, though the zero boundary is attracting 
and attainable, it is strongly reflecting — if the process reaches zero it leaves it immedi- 
ately and bounces back into the positive domain — see Revuz and Yor |48l p. 412]. We 
detailed in the introduction how this case is a major obstacle, particularly for direct 
discretization methods. A comprehensive account of direct discretization methods can 
be found in Lord, Koekkoek and Van Dijk 39 , to where the reader is referred. Based 
on our experience, the full truncation method proposed by Lord, Koekkoek and Van 
Dijk has so far proven to be the most accurate and efficient of this class. This is also 
evidenced by Andersen [5] and Haastrecht and Pelsser [22] who complete thorough 
comparisons. 

The method we propose follows the lead of Broadie and Kaya [TD] and Ander- 
sen [5], and is based on simulating the known transition probability density for the 
Cox-Ingersoll-Ross process. We quote the following form for this transition density, 
that can be found in Cox, Ingersoll and Ross [13], from a proposition in Andersen [5]. 

Proposition 1 Let F x 2^(z) be the cumulative distribution function for the non- 
central chi-squared distribution with v degrees of freedom and non-centrality parameter 
A: 



Set v :— Ak6/e 2 and define 



rj(h) 



4k exp(-Kh) 
l (l -exp(-Kh)) ' 



where h = t n +i — t n for distinct times f n +i > t n . Set A := Vt n ■ n{h). Then condi- 
tional on Vt n , Vt n+1 is distributed as exp(—Hih)/ri(h) times a non-central chi-squared 
distribution with v degrees of freedom and and non-centrality parameter A, i.e. 

P(V*„ +1 <x\V tn )= F xl{x) (x ■ v(h)/ exp(-Kh)). 

Andersen [S] and Haastrecht and Pelsser [55] produce very effective approximate 
sampling methods for the non-central X^(^) transition density. Our goal is to compare 
these methods with the Marsaglia generalized Gaussian (MAGG) method for sampling 
exactly from this transition density outlined in Section [5] Note that for Cox-Ingersoll- 
Ross sampling using the MAGG, we set A = Vt n ■ r/(h) and compute 

Vu+i = (n 2 + ■ • ■ + Yf N + Z 2 1 q + --- + Zl q ) ■ exp(-Kh)/r](h). 
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3.2 Asset price process simulation 

We follow the lead firstly of Broadie and Kaya [10], and then secondly, of Andersen [5]- 
By integrating the exact equation for the volatility from t n to t„+i an expression for 
f \/VrdWr is obtained. This is substituted into the exact equation for In St, itself in 
integral form between t n and t n +i- The result is the exact relation 

lnS tn+1 = \nSt n +j(V tn+1 -Vt n -K0h)+(e2—$) £ " + V r dr+y/l - p 2 ^ " + vT^dW* 

Since W 2 is independent of the process Vt we naturally have, in distribution, that 




where Z ~ N(0, 1). Now we make the only approximation, suggested by Andersen, to 
approximate the remaining integral in the expression for In St n+1 by the midpoint rule: 




Frdrw \h{V tn +V tn+X ). 



Using these last two replacements, and exponentiating, we arrive at the approximation 
Stn+i = exp(A + KxV tn + K 2 V tn+l + V K 3Vt n +K i V tn+1 - Z), 

where K = -hpn/e, K\ = h(np/e - 1/2) /2 - p/e, K 2 = h(Kp/e - 1/2) /2 + p/e and 
A3 = K4 = h{l — p 2 )/2. One caveat remains. The exact process St is a martingale, 
however using the presciption just given E[St n+1 \St n ] 7^ St„. However we can correct 
for this, we quote again from Andersen [5] p. 21]. 

Proposition 2 Let K\, K 2 , A3, K\ be defined as above. With s := K 2 + set 

M := m[exp(sVt n+1 )\V tn ] and A * := - InM - {K x + %K 3 )V tn . 

Then if we replace Kq by Kq in the scheme for St n+1 above, we have E[St„ +1 \St n ] ~ 
St n - 

Hence the final task is to compute M. Since we simulate Vt n+1 exactly we know 
M = E [exp(s-z)|Vt„] , where z ~ XvW, with v and A defined for the Heston model, and 
s = s-exp(— Kh)/n(h). Hence provided s < 5 we have M = exp(As/(l — 2s)) / (l— 2s) u ^ 2 . 
Consequently in our simulation for St n+1 , we take 

Kq = ~jz?28 + ^ ln(1 " 2§) ~ (Al + 2~ K3)Vt "- 

The requirement s < \ translates to a mild restriction on the stepsize h, which in 
practice is not a problem (see Andersen 5; p. 24]). 
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Case I 


Case II 


e 


1 


0.9 


K 


0.5 


0.3 


P 


-0.9 


-0.5 


T 


10 


15 


«(0), 


0.04 


0.04 



Table 1 Test cases from Andersen. In all cases u(0) = 100. 



3.3 Simulation results 

We test our Marsaglia generalized Gaussian method (MAGG) directly against Ander- 
sen's QE-M method. The performance of the method of Haastrecht and Pelsser [22] 
is similar to Andersen's; the reader interested in the actual comparisons is referred 
to their paper. We use two of Andersen's test cases for pricing long-dated European 
FX call options, maturing at time T with strike K; see Table [1] Let the exact option 
price at maturity be C. The error of the approximation is E — C — C, where C is the 
sample average of the simulated option payout at maturity. In our examples, we use a 
sample size of 10 6 . In Tabled we show the errors for the two test cases at three dif- 
ferent strikes K — 100, 140, 60; without any postprocessing such as variance reduction. 
In terms of accuracy the MAGG method competes very favourably with Andersen's 
QE-M method, as might be expected. In terms of efficiency, averaging over all stepsizes 
in case I, the ratio of total CPU times, of the MAGG method over Andersen's method 
was 0.4981, making MAGG two times faster. Averaging over all stepsizes in case II, 
the ratio is 0.8121. 



4 Concluding remarks 

We have seen that the method of sampling from the central xi distribution that we 
have introduced, based on sampling from the generalized Gaussian distribution, is 
more efficient than the acceptance-rejection method of Ahrens and Dieter. This is true 
for degrees of freedom v < 2 and in particular for values of v given to three or four 
significant figures. Combining this method for the central xt distribution, with a sample 
from the Poisson distribution with mean A/2, means that we can efficiently sample 
from a non-central xiW distribution, and therefore from the transition probability 
density of the Cox-Ingersoll-Ross process for the volatility in the Heston model. We 
demonstrated, that for values of the number of degrees of freedom v = 4«:6/e 2 in the 
Heston model, given to three significant figures, the method we proposed is in fact 
more efficient than the Andersen method for the same accuracy. An interesting issue 
that arises here is the structural stability, i.e. sensitivity to the calibrated parameter 
data, of the Heston model. For example, how sensitive are the computed option prices 
to having a fourth or fifth significant figure of accuracy? We intend to pursue this 
question next. 

In addition, Glasserman and Kim [20] have recently introduced a novel method for 
simulating the Heston model. As we can see from our analysis in Section [3T2l to compute 
the price process at the end-time T, we in essence need to sample from the distribution 
for J V T dr on the interval [0, T\. The transition density for this integral process over 
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Stcpsizc 


Case I 


Case II 


h 


Andersen Marsaglia 


Andersen Marsaglia 















Strike 


100 












1 


0.2211 


(0 


.012) 


-0.2374 


(0. 


.013) 


-0.4833 


(0. 


.042) 


-0.1404 


(0 


.042) 


1/2 


0.1164 


(0 


.013) 


-0.0707 


(0. 


.013) 


-0.0400* 


(0. 


.046) 


-0.0264* 


(0 


.044) 


1/4 


0.0143* 


(0 


.013) 


-0.0440 


(0 


.013) 


-0.0231* 


(0 


.044) 


0.0217* 


(0 


.048) 


1/8 


-0.0277* 


(0 


.013) 


-0.0050* 


(0 


.013) 


0.0807* 


(0 


.045) 


-0.0553* 


(0 


.052) 


1/16 


0.0162* 


(0 


.013) 


0.0019* 


(0 


.013) 


-0.0026* 


(0 


.042) 


0.0521* 


(0 


.046) 














Strike 


140 












1 


-0.0883 


(0 


.002) 


-0.0283 


(0 


.002) 


-0.3082 


(0. 


.036) 


-0.0926* 


(0 


.036) 


1/2 


-0.0274 


(0 


.003) 


-0.0121 


(0 


.002) 


0.0515* 


(0 


.040) 


0.0029* 


(0 


.037) 


1/4 


-0.0013 


(0 


.003) 


-0.0048 


(0 


.003) 


-0.0016* 


(0 


.038) 


0.0207* 


(0 


.043) 


1/8 


0.0047 


(0 


.003) 


-0.0011 


(0. 


.003) 


0.0740* 


(o. 


.039) 


-0.0327* 


(0 


.047) 


1/16 


0.0018 


(0 


.003) 


0.0015 


(0. 


.003) 


0.0069* 


(0 


.035) 


0.0509* 


(0 


.040) 














Strike 60 












1 


0.0317* 


(0 


.025) 


-0.1234 


(0 


.025) 


0.1180 


(0 


.048) 


-0.0379* 


(0. 


.049) 


1/2 


0.0345* 


(0 


.025) 


-0.0556* 


(0 


.025) 


0.1349 


(0 


.052) 


-0.0036* 


(0 


.050) 


1/4 


0.0111* 


(0. 


.025) 


-0.0388* 


(0 


.025) 


-0.0066* 


(0. 


.050) 


0.0290* 


(0 


.054) 


1/8 


0.0407* 


(0 


.025) 


0.0120* 


(0 


.025) 


0.0809* 


(0. 


.052) 


-0.0650* 


(0 


.058) 


1/16 


0.0284* 


(0 


.025) 


0.0003* 


(0 


.025) 


-0.0170* 


(0. 


.049) 


0.0492* 


(0 


.052) 



Table 2 Estimated error using 10 e paths. Sample standard deviations are shown in paren- 
thesis. As Andersen, we star results that are not statistically significant at the level of three 
sample standard deviations. 



the whole interval [0, T], given Vq and Vt, is well known and given in Pitman and 
Yor [47] . Its Laplace or Fourier transform has a closed form. Broadie and Kaya [10] 
use Fourier inversion techniques to sample from this transition density for J V T dr. 
Glasserman and Kim instead separate the Laplace transform of this transition density 
into constituent factors, each of which can be interpreted as the Laplace transforms of 
probability densities, samples of which can be generated by series of particular gamma 
random variables. The advantage of this method is that J V T dr is simulated directly on 
the interval [0, T\. However it does require approximation in the form of truncating the 
series, and the simulation of a sufficiently large number of gamma random variables. It 
would be of interest to compare the accuracy and efficiency of Glasserman and Kim's 
approach to ours, especially if we simulate the price process St directly as follows. 
We decompose J V T dr on [0,T] into subintervals [t n ,t n +x\, use a simple quadrature 
to approximate J V T dr on these subintervals much like Andersen, and simulate the 
transition densities required using the generalized Marsaglia method. We then only 
exponentiate at the final time T to generate an approximation for St (since we do not 
compute the price process at each timestep, this will be more efficient). However, one 
advantage of the approach we have taken in Section [372] for simulating the price process 
based on the method proposed by Andersen, is that it is more flexible. For example, 
it allows us to consider pricing path-dependent options. It would also be of interest to 
compare the accuracy and efficiency of this approach for the pricing of path-dependent 
options as well. 
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